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Abstract. We consider dynamically constrained Monte-Carlo dynamics and show that this leads to 
the generation of long ranged effective interactions. This allows us to construct a local algorithm for 
the simulation of charged systems without ever having to evaluate pair potentials or solve the Poisson 
equation. We discuss a simple implementation of a charged lattice gas as well as more elaborate off- 
lattice versions of the algorithm. There are analogies between our formulation of electrostatics and 
the bosonic Hubbard model in the phase approximation. Cluster methods developed for this model 
further improve the efficiency of the electrostatics algorithm. 
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1. Introduction 

Computer modeling of charged systems is demanding due to the range of the Coulomb 
interaction [1]. The direct evaluation of the Coulomb sum for N particles, Uc = 
J2i<j ^i^j I ^''^^ofiji requires computation of the separations between all pairs of par- 
ticles, which implies 0{N'^) operations are needed per sweep or time step. Most large 
scale codes in use at the moment solve Poisson's equation for the electrostatic potential 
using the fast Fourier transform [2] after interpolating charges to a grid. For an ensemble 
of charges interpolated to M nodes of a lattice this gives the electrostatic energy in a time 
which scales as 0{M\\iM). 

The use of a global solver for the Poisson equation leads to a strong bias in the choice 
of simulation techniques. Since the solution of the Poisson equation is unique the motion 
of a single charge requires the global re-calculation of the electrostatic potential, so that 
Monte-Carlo methods must (apparently) lead to hopelessly inefficient codes: All efficient 
large scale codes are based at the moment on molecular dynamics rather than Monte-Carlo. 
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2. Electrostatics can be formulated locally 

It is surprising that the Coulomb interaction poses such tremendous difficulty; after all 
the underlying Maxwell equations are local. The question then arises as to what part of 
Maxwell's equations give rise to electrostatic interactions. Is it possible to generate effec- 
tive Coulomb interactions in a manner better adapted to computer simulation? 

In electromagnetism Coulomb's law comes [3] from a local expression for the energy: 

C/ = |!2^d3r (1) 

where E is the electric field, and the imposition of Gauss' law 

div E - p/eo = 0. (2) 
This motivates the use of the following partition function for the electric field 



2({r J) = JveHS (div E - p ({r J) /eo) e'^/'^-^. 



(3) 



The charge density, p{r) — eiS{r — r^) and the charge of the i'th particle is e^. We 
change integration variables to e = E + V0p with V^0p = — p({ri})/eo and find 

2({ri}) = Jve <5(div e) e"^* /(-^^p)^ d^r^ 
We now notice that the cross term in the energy is zero by integration by parts: 

2U/eo = j {V(l)pf d\ + j e'^ d?v-2 j V(t>-e d\ (5) 

^ V ' 

=0 

so that 

Z = e-P'-^ n^M' J Ve S{dW e)e-^* (6) 

We conclude that Z{{ri}) = Zcouiomb{{^i}) x const: Relative statistical weights are 
unchanged if we allow the transverse field to vary freely, rather than quenching it to zero. 
To sample the partition function Eq. (3) we only have to find solutions to Gauss' law. If 
we dispense with the condition curl E = usually imposed in electrostatics, E can be any 
solution to the equation div E — p/eo = 0; the general solution is E = — grad + curl Q 
with Q arbitrary. 



2.1 Monte-Carlo sampling 

We have formulated a Monte-Carlo algorithm [4] that uses the Metropolis method together 
with the energy Eq. (1). In order to generate configurations according to the statistical 
weight of Eq. (3) we need to 

• move particles without violating the constraint of Gauss' law to preserve the delta- 
function on div E — p/eo; 

• integrate over the transverse degrees of freedom e, of the electric field. 
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Local Simulation Algorithms 
Discretization 

The system is discretized by placing charged particles on the M = vertices of a periodic 
cubic lattice, {i}. The components of the electric field Ei j are associated with the 3Af 
links {i,j} of the lattice. There are 3 A/ plaquettes on the lattice each defined by four 
links. We use the notation Ei^2 to denote a local contribution to electric flux leaving 1 
towards 2. If the link from 1 to 2 is in the positive x direction we consider that this is the 
local value of the x component of the field E. The 3A/ variables Eij are thus grouped into 
M three dimensional vectors. 



Particle motion 



l/4ii 



e=l 



1/4 



l/4v 



1/4 



1/12V 



1/12 



Figure 1. A local solution to Gauss' law. Numbers on each link correspond to the 
electric flux in the direction of the arrows. There is a charge e/eo ~ 1 on the leftmost 
node of the figure. 



Two nodes of the lattice are shown in Fig. (1). Gauss' law is interpreted in the integral 
form / E • dS = e/eo where e is the charge enclosed by a surface enclosing each site. On 
our lattice we associate the fluxes with the links; if a charge of e/eo = 1 is at the leftmost 
site we can satisfy Gauss' law with a flux of 1/4 on each outward link. On the rightmost 
node the charge is zero and the signed sum of the fluxes at this site is also zero. 
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Figure 2. After moving the particle from the leftmost to the rightmost node we wish 
to find a new solution to the Gauss constraint. This is done by modifying the flux on a 
single link connecting the two sites from 1/4 to —3/4. 
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We now displace the particle to the rightmost site: Fig. (2). We must find a new solution 
to Gauss' law. One solution is found by taking the field on the link between the sites, E and 
modifying it with the law E ^ E — e/eo ot from 1/4 to —3/4. This is the essence of the 
method: Gauss' law allows an entirely local modification of the fields, whereas Poisson's 
equation implies a global update of the potential at all sites of the lattice. The energy 
change, AU ~ eo ((—3/4)^ — (1/4)^) /2, associated with the displacement is used in a 
standard Metropolis algorithm to decide whether to accept the update. 



Transverse field 

The update of the field, slaved to the particle motion is not in itself sufficient to fully sample 
the partition function, Eq. (3); we must also integrate over the transverse field e. One can 
group four links into plaquettes and modify the four Unks by the same increment, A so as 
to conserve Gauss' law: Fig. (3). In our code particle motion and plaquette updates are 
mixed stochastically with probabilities p and (1 — p). We chose A uniformly distributed 
in [— Ao,Ao], where Aq is chosen so that the Monte-Carlo acceptance probability of this 
move is close to 1/2. 



-4,3 



'3,2 
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Figure 3. The circulation of a plaquette is incremented by 4A by modifying the field 
ofeachlink: £4,3 £;4,3-A, £4,1 £;4,i-|-A, £3,2 ->■ £3,2-A, Si,2 ->■ Ei^2+A 



2.2 Large scale dynamics 

The above algorithm is based on local, stochastic updates for both the particle positions 
and the electric field. It is thus rather natural that both the particles and the electric field 
exhibit diffusive dynamics. It can be shown [6] that the field dynamics are governed by a 
Langevin equation: 

— = (V^E - grad p/eo) - J/eo + curl ({t, r) (7) 

which replace Maxwell's equations in our Monte-Carlo scheme. The diffusion coefficient 
of the transverse electric field is related to p and Aq. The analytic structure of Eq. (7) 
is particularly interesting when there are free, mobile charges so that the electric current is 
linked to the electric field via the equation J = ctE. In this case the dispersion law of the 
electric field develops a gap so that 
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Local Simulation Algorithms 

iwqAE= (CT/eo + -D£;q^)qAE. (8) 

The relative diffusion rates of the electric and charge degrees of freedom are adjusted by 
modifying p. We initially believed that the algorithm should be run in a regime in which 
the effective diffusion coefficient of the electric field is comparable to or somewhat larger 
than that of the particles: The particles are then always interacting via a field that is close to 
equilibrium and are not slowed down by the "drag" from the electric degrees of freedom. In 
dilute systems this choice would mean that much more CPU time must be spent updating 
plaquette degrees of freedom than those corresponding to the particles and the algorithm 
would become inefficient. 

Simulations show that one can use a much lower rate of update for the plaquettes. There 
are several ways of understanding the efficiency of the code [7]: 

• The gap in the electric dispersion law Eq. (8) gives fast relaxation of the electric 
degrees of freedom even at long wavelengths and when De is small. 

• Motion of particles akeady integrates over the electric field: Motion of a particle 
around a single plaquette has exactly the same effect as an update of a plaquette by 

A = e/eo. 




• Particles only, T=0. 15 
■ Particles and plaquettes, 
O Particles only, T=0.5 
Q Particles and plaquettes, 
y = x/(l+x) 



(1-kV/8)C0^/k^ 



20 



Figure 4. Charge structure factor as a function of ujq. At high temperatures the same 
structure factor Eq. (9) is found with and without plaquette updates. At lower tem- 
peratures this law is still observed when the plaquettes are updated; without plaquette 
updates one finds S ~ Ljq. Fit to d = 1.29. L — 15, N — 336. 



We simulated the above lattice gas and measured the charge-charge structure factor. 
Debye-Hiickel theory applied to a lattice gas gives the following expression for S{q) = 
ijjql[K} + tjq) with ujq = Y^i=\ ^(1 ~ cosqi) and — ce^ /ksTeQ, with c the charge 
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density. When q ^ then ujq q^. When the finite size of particles is also taken into 
account [8] the long wavelength structure factor must be slightly modified. We thus fit our 
numerical structure factors to the expression 

^^^^ = K^+u:,i^-KH^/B,) 

where d is the particle diameter. In Fig. (4) we show results of simulations at kBTeo/e^ = 
0.5 and kBTeo/e^ = 0.15 which are in very good agreement with this expression. 

Numerical experimentation leads to a rather remarkable result: If one drops the plaquette 
updates entirely and works at a sufficiently high temperature the charge-charge correlation 
functions are almost indistinguishable, Fig. (4), from those in which the system is sim- 
ulated with plaquette updates. The system is clearly not fully ergodic, since only values 
of the plaquette circulation which are integer multiples of e/eo are possible. Yet charge 
correlation functions appear to be very similar to those of a fully equilibrated system. 

It is only at lower temperatures that there is a qualitative change in behaviour and the 
system has very different behaviors depending on whether the plaquettes are independently 
excited or not. One possible explanation is that without plaquette updates the charges are 
condensing into tightly bound, neutral pairs. Indeed, if one studies a system with just two 
charges at low temperatures one finds a finite, constant line tension between charges, and 
thus a potential between two charges which increases linearly with separation. 



3. Improved discretization 

The lattice gas described above is not yet useful for simulations in condensed matter 
physics: The algorithm becomes very inefficient at low temperatures. Motion of a par- 
ticle between two nodes leads to a finite modification of the field on the connecting link, 
and thus a finite energy barrier; the Monte-Carlo acceptance rate falls off exponentially 
at low temperatures. In practice this leads to hopelessly inefficient simulations for tem- 
peratures such that /csTeo/e^ < 0.15. In condensed matter physics it is usual to work 
in terms of a normalized Bjerrum length, £b- the length at which the electrostatic inter- 
action between two particles is equal to fc^T when measured in units of particle size: 
ts/a = e'^/AncokBT. Typically we are interested in ~ 5 - 20. /csTeo/e^ = 0.15 
corresponds to only £b ~ 1/2, we need to gain at least a factor of ten in temperature to 
have a useful code. 



Charge spreading 

One can avoid the fall off in efficiency at low temperatures by spreading the charge over 
several lattice sites. In Fig. (5) we have plotted the Monte-Carlo acceptance rate for particle 
motion as a function of temperature for a system of dimensions L = 15 with two free 
charges. We associate the charge of each particle to sites with n = 1,2,3. We find 
that the temperature at which the code becomes inefficient varies as T ^ When we 

measure the Bjerrum length in terms of the size of the spread-out particle we find that we 
can simulate at a Bjerrum length £b < without the acceptance rate falUng too low. 
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Figure 5. Evolution of acceptance rate as a function of the degree of spreading. 
Dashed lines: guide to eye. Temperature in units of e^/eofes. 



Off-lattice methods 

By allowing the particles to move in the continuum and then interpolating the charges to 
the network (as is usual in Fourier codes [2]) we can continuously tune the Monte-Carlo 
step size so that the acceptance rate does not fall off with temperature. Many choices 
of interpolation scheme are possible, in our codes [5] we have used a scheme based on 
piecewise quadratic splines, leading to interpolation to a cube of 27 sites around each 
particle. 

When a particle moves in the continuum we need to generalize the update rule E — > 
E — e/eo that generates the local updates to the electric field. We must generate an electric 
current J which is a solution to the continuity equation div J — —Ap, where Ap is the 
finite variation of the charges on the nodes calculated from the spline interpolation scheme. 
We then update the electric field with E — *■ E — J/eg. The equation for J is clearly under- 
determined. We use this fact to our advantage to find one solution which is simple. This is 
illustrated in Fig. (6). We construct a Hamiltonian path which visits each node where the 
interpolated charge has changed. The current is then calculated on each link as the sum of 
the charge variations on all previous nodes on the path. 



4. Cluster sampling for the electric field 

In many physical situations the algorithm as described above is remarkably efficient: The 
presence of a finite concentration of free charges always gives rise to a finite a and thus 
fast field relaxation: the slowest modes in the system are then modes of density relaxation. 
However, one is sometimes interested in highly heterogeneous systems where charges are 
either not mobile (perhaps they are attached to a polymer or a surface) or are excluded from 
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Acj ; Ae^ \ Ae3 

J=-Aej/8Q J=-(Aej+ Ae2)/eQ 
Figure 6. Illustration of charge interpolation and the current update using a Hamil- 
tonian path in two dimensions. A charge (open circle) interpolates onto the 3^ = 9 
lattice sites (solid). The partial charges on these sites change by an amount Aej when 
the particle moves. The current J flows along the Hamiltonian path and terminates at 
the last site, since Aej = 0. 



certain parts of the simulation volume (perhaps due to the presence of a large colloidal 
inclusion). Similarly in dipolar fluids such as water in the absence of free charges ct = 0. 
In this case we have to rely on the slower diffusive motion of the electric field in Eq. (7) 
driven by the plaquette updates. We will now show how to generate updates which give 
rise to fast field equilibration even for the case cr = using a cluster algorithm for the 
electric field variables. 

There is a remarkable analogy between the constrained partition function Eq. (3) and the 
bosonic 2+1 dimensional Hubbard model at zero temperature in the phase approximation 
L9J. Its partition function is written in terms of a constrained vector field discretized on a 
lattice: 

E ^^P(-^ / ^'^'0 n'^(divE). (10) 

E integer \ J / r 

The major difference compared with the partition function Eq. (3) is that the field variables 
on each Unk are now integers rather than reals; the field variables represent the real bosonic 
currents in the system. K is the ratio between kinetic energy and the repulsive energy and is 
interpreted as an effective temperature. The partition function Eq. (10) has been extensively 
studied [9] and has a phase diagram with a small K (large repulsion) insulating phase and 
a superfluid phase at large K (large kinetic energy) .There is a phase transition between 
these phases at K = 0.33305(5). We have already seen an example of a similar partition 
function when we suppressed the plaquette updates in the electrostatic algorithm. The very 
different structure factors observed in Fig. (4) are presumably Unked to the two phases of 
the Hubbard model. 

Early numerical studies of Eq. (10) were based on local updates to plaquette degrees of 
freedom, in a manner very similar to that described above for the electrostatic algorithm. 
More recently much more efficient algorithms have been found based on cluster or worm 
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sampling [10]. The algorithm works by nucleating a pair of positive and negative pseudo- 
particles on the same site of the system. One of the pseudo-particle is then moved on the 
lattice so that div E 7^ 0. The motion of this pseudo-particle is a random walk biased by the 
energy in Eq. (10) [10]. The pseudo-particle eventually returns to its stationary partner and 
the two particles then annihilate each other. When the pseudo-particle moves, the electric 
field is updated exactly as described above for the electrostatic algorithm, E ^ E — e/eo. 
At the moment of annihilation the total weight of the generated path or cluster is compared 
with its time reversed version and the update is either globally accepted or refused [10]. 
The algorithm is very reminiscent of the creation of virtual pairs of particles in quantum 
mechanics. In the original version of the algorithm [10] the pair of pseudo-particles had a 
charge e/eo — ±1 in order to generate fields with integer flux on the links. 
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Figure 7. Decay rate, F of the transverse electric field as a function of ujq. □: empty 
box, plaquette updates only illustrating diffusive field motion F — De^^q', time scale: 
L'' Monte-Carlo tries. ^1^: two particles in box with plaquette updates F = a /e+Dsi-L'q', 
time scale: Monte-Carlo trials, 1:1 split in Monte-Carlo attempts between particles 
and plaquettes. O: cluster algorithm showing uniform relaxation across all modes; time 
scale: a single cluster update. L = 15, fc_gTeo/e^ — 0.5 . 

Given the efficiency of this algorithm in sampling Eq. (10) we have explored the interest 
of the method for sampling the electrostatic partition function Eq. (3). In the electrostatic 
problem we are interested in integrating over all values of the electric field. Thus we 
create the pair of virtual particles with a random charge distributed according to a uniform 
distribution between — e,„ and Cm- If is small we find that the typical cluster is very 
large containing about 0.7 x i"^ links. In order to maximize the rate of change in the field 
on the links we should use a large value for £„. However, the use of a too large value 
for Cm prevents the propagation of the virtual particle, which backtracks too readily in its 
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own path and only visits 0(1) sites before the annihilation. This is again understood in 
the boson analogy as being due to the different behaviour of a particle in a superfluid or 
insulating phase. In practice we use ^ fesTeo- 

In Fig. (7) we plot the decay rate for the Fourier components of the transverse electric 
field for three different simulations. The first simulation is a box without free charges sim- 
ulated using simple plaquette updates. As predicted by Eq. (7) the dynamics of the field is 
diffusive, T = Dsq^ iov small q. We then add two free charges to the system and measure 
the same relaxation rate. We see that a finite gap is introduced into the dispersion law 
as predicted by Eq. (8). Finally we simulate the system without charges using the cluster 
algorithm. We find that the relaxation rate is weakly dependent on the mode and we have a 
method which equilibrates the transverse electric field in about two cluster moves per Unk, 
rather than the O(L^) sweeps needed for the diffusion of the electric field with local pla- 
quette updates. Mixing particle motion with occasional cluster updates in a heterogeneous 
system now allows one to ensure fast relaxation of the electric degrees of freedom even in 
situations where a is zero. 

We understand the efficiency of the cluster algorithm in both the Hubbard model and 
in the electrostatic algorithm as being due to the dispersion law Eq. (7) in the presence 
of a current J = ctE. We have seen that if a is non-zero then we develop a gap in the 
electric spectrum and the system is described by a dynamic exponent z = Q rather than 
z = 2 characteristic of diffusion. In the cluster algorithm the finite conductivity comes 
from the virtual particles rather than from true free ions but their effect is very similar on 
the relaxation spectrum for the electric field. 
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